Spin-dynamics simulations of the triangular antiferromagnetic XY model 



(N 
O 

o 

(N 



J2 

o 

CD 



Kwangsik Nho and D. P. Landau 

Center for Simulational Physics, University of Georgia, Athens, Georgia 30602 

(Dated: February 1, 2008) 

Using Monte Carlo and spin-dynamics methods, we have investigated the dynamic behavior of the 
classical, antiferromagnetic XY model on a triangular lattice with linear sizes L < 300. The tem- 
poral evolutions of spin configurations were obtained by solving numerically the coupled equations 
of motion for each spin using fourth-order Suzuki- Trotter decompositions of exponential operators. 
From space- and time-displaced spin-spin correlation functions and their space-time Fourier trans- 
forms we obtained the dynamic structure factor 5(q, w) for momentum q and frequency u). Below 
Tkt (Kosterlitz-Thouless transition), both the in-plane (S xx ) and the out-of-plane (S zz ) components 
of S(q,u>) exhibit very strong and sharp spin- wave peaks. Well above Tkt, S xx and S zz apparently 
display a central peak, and spin- wave signatures are still seen in 5 122 . In addition, we also observed 
an almost dispersionless domain- wall peak at high uj below T c (Ising transition), where long-range 
order appears in the staggered chirality. Above T c , the domain-wall peak disappears for all q. The 
lineshape of these peaks is captured reasonably well by a Lorentzian form. Using a dynamic finite- 
size scaling theory, we determined the dynamic critical exponent z — 1.002(3). We found that our 
results demonstrate the consistency of the dynamic finite-size scaling theory for the characteristic 
frequeny uj m and the dynamic structure factor S(q,w) itself. 
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I. INTRODUCTION 

The study of the dynamic properties of classical spin 
systems has been extensively carried out using both theo- 
retical and simulational approaches. Irutheir work on the 
theory of dynamic critical phenomenal^, Hohenberg and 
Halperin proposed a number of different dynamic uni- 
versality classes based upon the conservation laws. The 
dynamic critical behavior is describable in terms of a dy- 
namic critical exponent z, which gives rise to different 
dynamic universality classes: z depends on the conser- 
vation laws, lattice dimension, and the static critical ex- 
ponents. Now spin dynamics simulations have become a 
mature method for probing the time dependent behav- 
ior of magnetic systems (see Ref. 2 for a recent review). 
This approach has resulted in high-quality dynamic criti- 
cal exponent estimates for numerous models. In a recent 
high-resolution spin-dynamics study,B Tsai et al. made 
direct, quantitative comparison of both the dispersion 
curve and the lineshapes obtained from their simulation 
data with recent experimental results for RbMnF3 us- 
ing the Heisenberg antiferromagnet model and a newly 
developed approach|-pf higher order decomposition time 
integration methodsa and obtained a good agreement. 

In this paper we consider the classical two-dimensional 
antiferromagnetic XY model on a triangular lattice 
(TAFXY). This frustrated 2D spin syste 
ceived much attention during the last decadeE 
The TAFXY model displays rich low-temperature phase 
structures and critical phenomena because frustration in- 
troduces additional discrete symmetries resulting from 
the chiral degrees of freedom. The ground states of this 
model are composed of three interpenetrating sublattices 
with lattice vectors of length v3- Spins on each sub- 
lattice are ferromagnetically ordered and spins on dif- 
ferent sublattices are oriented ±27r/3 with respect to 
each other. The TAFXY has a continuous U(l) sym- 



metry associated with global spin rotations and a dis- 
crete Z2 symmetry because of the double degeneracy 
of the ground-state chirality configurations. This model 
has two order parameters: the staggered in-plane mag- 
netization and the staggered chirality. There have been 
ongoing controversies concerning phase transitions and 
the nature of the phase transitions. SomeEffil suggested 
that the system undergoes a single phase transition, 
generalW|-|i4eMing non-Ising critical exponent; however, 
othersQt1all3'yj supported a double phase transition sce- 
nario with a Kosterlitz-Thouless(KT)-like transition fol- 
lowed by a second-order chirality-lattice melting transi- 
tion at a slightly higher temperature. In addition to spin 
waves and vortices, which are fundamental excitations 
in the ferromagnetic or antiferromagnetic XY model on 
a bipartite lattice, the triangular antiferromagnetic XY 
model has another excitations associated with domain- 
wall formation between two different ground states in 
view of the chirality configuration of ground statesa. 

A good experimental realization of the frustrated an- 
tiferromagnetic XY model on a stacked triangular lat- 
tice is ABX9 type halides (like CsNiCi3, CsMnBr3, 
CsCuCl3)E3lij. The magnetic B 2+ ions form a triangular 
antiferromagnet within the ab planes and a ferromagnetic 
coupling along the c direction produces three dimensional 
order. 

While the static critical behavior of the TAFXY model 
has been the subject of much interest during the last 
decade, the theoretical study of dynamic critical behavior 
has not been carried out for the TAFXY model. 

In the present work, we have studied the dynamic be- 
havior of the TAFXY model using Monte Carlo and spin- 
dynamic methods(MC-SD). Although in the usual planar 
XY model each spin has two components, here we use 
the three-component XY model in order to study the real 
dynamics. These-4w© models belong to the same static 
universality classtilta, but the planar XY model will not 
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have true dynamical behavior. Ill 

The estimates for Tkt and T c used here are taken from 
a recent high-precision MC study where the critical tem- 
perature T c associated with the chirality phase transition 
and the KT transition temperature Tkt associated with 
unbinding of vortex pairs have been determined as T c — 
0.412(5)J/fe s and T KT = 0.403(1) J/k B M 



II. MODEL AND SIMULATION METHOD 

We will first briefly describe the model and the numer- 
ical method used here and show how the dynamic struc- 
ture factor is computed. The classical antiferromagnetic 
XY model is described by the following form: 



S tt (q,w) = ye ,(1 '( r - r ') / e iujt C kk (r-r',t) — , (3) 

where k = x, y, or z is the spin component and the space- 
displaced, time-displaced spin-spin correlation function is 
given by 

C kk {r-v',t-t') = (S k (t)S k ,(t')) 

-(S k (t))(S k (t<)). (4) 

The dynamic critical exponent z can be determined by 
using the dynamicJimte-size scaling theory developed by 
Chen and LandauEj: 



<tj> 



(1) 



where the summation is over all nearest neighbors, Sj 
= (S?,S?,S*) is a three-dimensional classical spin of 
unit length at site i, and J is the positive (antiferromag- 
netic) coupling constant between nearest-neighbor pairs 
of spins. We consider triangular lattices of size L x L 
along the primitive vector directions^, ^e x + ^i y ), 
containing N = L 2 sites and 27V elementary triangles. 
In our calculations, periodic boundary conditions are ap- 
plied along the primitive vector directions. 

We use a hybrid Monte Carlo procedure which, con- 
sists of a combination of the—Metropolis updateO and 
the over-relaxation algorithnxZl ( see also Ref. 18). One 
hybrid Monte Carlo step consists, of two Metropolis and 
eight over-relaxation updateal3'E2l. Using this hybrid al- 
gorithm, we generated approximately 700-8500 equilib- 
rium configurations at a given temperature. Typically 
1000 hybrid Monte Carlo steps were used to generate 
each equilibrium configuration, which is then evolved us- 
ing the equations of motion given by: 



xf(q) 



G kk (u,L z ,qL), 



(5) 



where we do not introduce a frequency resolution func- 
tion in order to smoothen the effects of finite length of 
time integration because of very long integration times 
used in our spin dynamics. Xl'°( ( i) is the total integrated 
intensity given by 
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(6) 



The characteristic frequency ^^(qL) is a median fre- 
quency determined by the constraint: 



Xf(q) 



(7) 



In the dynamic scaling theory, the finite-size scaling 
expression for the median frequency uj kk is given by 



Starting from a particular initial spin configuration, we 
performed numerical integration of these equations of 
motion using a recently developed 4th-order Suzuki- 
Trotter decomposition method], and the integration is 
carried out to a maximum time t max , typically of the 
order of t ma xJ = 680 and t cu toffJ = 600, with a time 
step of St = 0.2J -1 . We compute the thermal average 
of a time-dependent observable by averaging over all the 
values of the observable obtained by evolving all the in- 
dependent initial equilibrium configurations. 

The dynamic structure factor S kk (q,u>) is the space- 
time Fourier transform of the position and time displaced 
spin-spin correlation function. It is defined for momen- 
tum transfer q and frequency transfer u> as followaiij: 



uj kk = L- z n(qL). (8) 

Using Eq. (@), we estimate the dynamic critical exponent 
z from the slope of a graph of log(w^) versus log(L) at 
fixed value of qL, and we test the dynamic finite-size 
scaling theory using Eq. (||). 

III. RESULTS AND DISCUSSION 

We first compare the behavior of the dynamic structure 
factor S kk (q,uj) at several different temperatures, where 
k refers to the x (in-plane) or z (out-of-plane) component 
and we were limited to the [11] reciprocal lattice direc- 
tion, i.e. q = (q, q) with q determined by the periodic 
boundary conditions, 
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FIG. 1: Temperature dependence of the dynamic structure 
factor S kk (q,io). L=240 and q= 2tt/15 in the [11] direction 
in all cases: (a) in-plane component; (b) out-of-plane compo- 
nent. 
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Figure 1 shows the temperature dependence of the 
dynamic structure factor S kk (q,ui) for a lattice size of 
L=240 and a particular momentum g=27r/15. For T < 
Tkt our results for the dynamic structure factor S xx 
show a very strong spin-wave and a central peak which 
are visible as pronounced peaks at the spin-wave fre- 
quency u>(q) and at u> = 0, respectively. .S 122 has struc- 
tures with less intensity than S xx (note the change in 
scale between Figs. 1(a) and 1(b)). There is a sharp 
spin- wave peak and no central peaks are visible in S zz . 
The positions of the spin- wave peaks are the same for S xx 
and S zz and as the temperature increases, the spin-wave 
peak broadens, its position moves towards lower u>, and 



FIG. 2: Temperature dependence of the dynamic structure 
factor S xx (q,ui) at high frequencies. L=240 and q= 27r/15 in 
the [11] direction in all cases. For clarity, a frequency resolu- 
tion function is used to smoothen the oscillations due to the 
finite t max . We can see a domain- wall peak for T < T c . The 
peak broadens as the temperature increases. However, the 
domain- wall peak disappears at T > T c . 
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FIG. 3: Lattice size dependence of dynamic structure factor 
S xx (q,uj) around the domain- wall peak position at fixed mo- 
mentum q— 2-7T/15 in the [11] direction at T = 0.5Tkt- The 
intensity of the domain-wall peak does not depend on lattice 
size and the data essentially collapse onto a single curve. 



its amplitude decreases. Above Tkt, S xx and S zz ap- 
parently display a central peak, but spin-wave signatures 
are still visible in S zz . The same qualitative behavior was 
observed in spin-dynamics simphitions of the 2D classical 
XY model on a square latticeEfl. n 

As has been pointed out by Lee et alB, in addition 
to spin waves, there is another type of elementary ex- 
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citation associated with domain-wall formation between 
regions with opposite staggered chirality. In order to ob- 
serve the dynamic effects of the destruction of chirality 
order, we investigate the temperature dependence of the 
dynamic structure factor S xx at high frequencies. We il- 
lustrate this in Fig. 2. The intensity of the domain- wall 
peak is 10 -2 of that of the spin-wave peak at small q ( 
note the change in scale between Figs. 1(a) and 2.). The 
relative intensity decreases at large q. For clarity, a fre- 
quency resolution function with Su>= 1.2ir/t cu t ff is used 
to smoothen the oscillations due to the finite t max . Below 
T c , where long-range order appears in the staggered chi- 
rality, we observe an almost dispcrsionless domain-wall 
peak at high ui typical of an Ising model. As the temper- 
ature approaches the chirality transition T c from below, 
the domain-wall peak broadens and its relative intensity 
decreases. The domain-wall peak disappears for all q 
above T c . This result indicates the presence of an Ising- 
like phase transition connected with the loss of chirality 
order. Fig. 3 shows S xx (q,w) for three different system 
sizes with the same value of q — 27r/15 at high lo in or- 
der to see the finite-size effects of the domain-wall peak. 
The intensity of the domain-wall peak does not depend 
on lattice size, the data essentially collapse onto a sin- 
gle curve, and finite-size effects are not visible, whereas 
the intensity of the spin-wave peak depends strongly on 
lattice size. 

The in-plane dynamic structure factor S xx (q,Lo) for 
the 2D XY model below the topologicaLtransition tem- 
peratupj Tkt was analyzed by yillainEJ, Moussa and 
VillainE 2 ], and Nelson and FisherEil They found that 
S xx (q, lo) has spin- wave peaks diverging at the spin- wave 
frequency Lo p (q). Menezes et alE3 calculated 5'(q, lo) and 
found a diverging spin-wave peak similar to that of Nel- 
son and Fisher and a logarithmically, diverging central 
peak. However Evertz and Landaui-3 found that their 
MC-SD data for the shape of S(q, lo) is not well described 
around the spin-wave peak by the above theoretical pre- 
dictions, rnm 

Mertens et aZ.EjO calculated S(q : lo) above Tkt and 
found a Lorentzian central peak for S xx and a Gaussian 
central peak for S zz . Very recently, Wysin et alfcjl cal- 
culated S*(q, lo) and found nondivergent spin-wave peaks 
and weak peaks on both sides of the spin-wave peaks. 
Below Tkt, these theoretical approaches motivated us 
to fit the lineshape of the dynamic structure factor S xx 
to a Lorentzian form. We found that the lineshape of 
S xx is reasonably well captured by a multiple Lorentzian 
form: 
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where the first term corresponds to the central peak, the 
next two terms are from the spin-wave creation and an- 
nihilation peaks at lo=±lo s , and the last two terms are 
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FIG. 4: In-plane component S xx of the dynamic structure fac- 
tor for L = 300. (a) T= 0.9T KT ; (b) T= Tkt at q=2-K /75 in 
the [11] direction in all cases. The symbols represent simula- 
tion data and the solid line is a fit with the multiple Lorentzian 
function given in Eq. (lid). 



contributions from the domain-wall peaks at u—±ujd. In 
order to fit the simulated lineshapes to the Lorentzian 
form Eq. (]Io|), we used a frequency range < lo < 5. 
We find that the Lorentzian lineshapes fit our simulation 
data reasonably well. Illustrations of the fits using Eq. 
( pj]| ) to the simulated lineshapes at T— O.QTkt and at 
Tkt are shown in Fig. 4. In the figure the symbols rep- 
resent our simulation data and the solid line is a fit with 
the Lorentzian function given in Eq. (|Io|). In Fig. 4(a) 
and (b), we find additional small peaks on both sides 
of the spin-wave peak. One simple explanation which 
is consistent with the data is the presence of two-spin- 
wave peaks. We could see that the position of these ex- 
tra peaks corresponded to frequencies of two spin-wave 
addition and difference peaks (marked as (1), (2), and 
(3) in Fig. 4(a)). These dynamical features have been 
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FIG. 5: Dispersion relation uj p (q) of the spin-wave peak in 
S xx (q,co) for the [11] direction as a function of momentum. 
The inset shows the dispersion relation w<j(g) of the domain- 
wall peak. 
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FIG. 7: Finite-size scaling of the median frequency( ujf L z as 
a function of qL z ) at (a) T = 0.9T KT ; (b) Tkt using 2=1.0. 
The data for u>^ show good scaling behavior for T < Tkt- 
Statistical errors of the data are smaller than the symbol sizes. 



FIG. 6: Finite-size scaling plot for the median frequency ai^f 
at T = Tkt- The dynamic critical exponent value obtained 
from the slope of the log-log plot is z—1. 000(1) for n=l using 
all lattice sizes and 2=1. 003(2) for n—2 using the largest four 
lattice sizes. The solid lines display a fit to a straight line. 
Statistical errors of the data are smaller than the symbol sizes. 



noticed in earliex_MC-SD simulationsHii-ffia and a theo- 
retical approacr£3. As we did not generalize Eq. (|i~0| ) to 
include the extra peaks, this extra structure causes the 
lincshape to depart from a Lorentzian form. 

The dispersion relations u) p (q) of the spin- wave peak 
and u>d(q) of the domain- wall peak can be obtained 
from the above Lorentzian fit. Fig. 5 shows the dis- 
persion curves as a function of momentum. The inset 
shows the dispersion relation u>d(q) at T = 0.5Tkt and 



T = Q.9Tkt, where obtaining a good u>d(q) is difficult 
because of the small peak height and the broadening of 
the domain- wall peak (see Fig. 2). As expected, u) p (q) is 
linear for small q and u>d(q) is approximately constant. 

The dynamic critical exponent z can be extracted from 
the finite-size scaling of the characteristic frequency uj mi 
as described in a previous section. We calculated uj^(qL) 
using S xx with qL=4ir(n=l) and qL=8ir(n=2) for several 
lattice sizes at Tkt- The log-log plot of the characteristic 
frequency w^f as a function of the lattice size L is shown 
in Fig. 6 where the estimated error bars for individual 
points are smaller than the symbol sizes. The solid lines 
display a fit to linear behaviors. From the slope of a 
log(w^f ) vs log(L) graph, we obtained the dynamic crit- 
ical exponent z — 1.000±0.001 for n—1 using all lattice 
sizes and z= 1.003±0.002 for n—2 using only the largest 
four lattice sizes. Within their respective error bars, the 
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two estimates for z agree. Our final estimate for the dy- 
namic critical exponent is z = 1.002(3), where the error 
bar reflects the fluctuations in the different estimates of z. 
For comparison, we proceeded to estimate the dynamic 
critical exponent z at T c . We obtained the dynamic crit- 
ical exponent z=1.051±0.005 for n=l (not shown). We 
then calculated ui kk (qL) using S kk for all q in the [11] 
direction and all lattice sizes L and we graphed lo^L z 
as a function of qL z for T < Tkt using z=1.0. This 
graph is shown in Fig. 7. In the figure, statistical errors 
of the data are smaller than the symbol sizes. As we 
know, the dispersion curve flattens for each finite lattice 
size when q becomes large; therefore, the data start to 
move away from the asymptotic behavior at progressively 
larger values of qL. The data for ui^ show good scaling 
behavior for T < Tkt- The out-of-plane characteris- 
tic frequency lu^ has the same scaling behavior as the 
in-plane component at temperatures T < Tkt- Interest- 
ingly, at T = T c , oj^ also shows the same good scaling 
behavior with z=1.0; however, we do not observe similar 
scaling behavior in w^f at T = T c . The same behavior 
was observed by spin-dynamics simulation of the-2-D XY 
model on a square lattice in a qualitative sense.EJ 

We had also attempted to estimate the q dependence 
of the half- width r fefc (q) of spin- wave peak in the critical 
region. Here we expect the simple form T kk (q) ~ q z . The 
log-log plot of the half- width r fefc (q) from our simulations 
is shown in Fig. 8. For T zz (q), we estimated the error in 
the fitted parameters by fitting the lineshapes using three 
different ranges of frequency around the spin-wave peak. 
The fitted parameters varied when different frequency 
ranges were used in the fit. For F xx (q), however, as there 
are three peaks( central, spin-wave, and domain-wall ) 
in S xx (q,Lo) below T c , obtaining a good fitted parameter 
was more difficult than T zz (q). First, the lineshape pa- 
rameters are obtained from the fit to only one frequency 
range. Using these parameters, the central peak and the 
domain-wall peak are subtracted. The error is estimated 
by varying the size of the frequency window around the 
spin-wave peak. In contrast to our expectation, we ob- 
served two different dynamic exponents for two different 
regions of q in T kk (q) at T — 0.9Tkt- In order to check 
our result, we calculated T zz (q) for the two-dimensional 
XY model on a square lattice and obtained two different 
exponents as well(not shown). 

In order to test the dynamic finite-size scaling theory 
of the dynamic structure factor S kk itself through the 
use of Eq. (H), we plotted S kk (q,uj)/L z x'l k (q) vs uiL z 
for four lattice sizes with the dynamic critical exponent 
set to z=1.0. For gL=47r (n=l) and T—Tkt, the result- 
ing scaling plot for S xx (q, to) is shown in Fig. 9. Data 
points collapse onto the same curve for sufficiently large 
L, namely lattice sizes L > 240. Only the data from 
small lattices deviate systematically. Scaling works quite 
well for all frequencies for lattice sizes L=240 and 300. 
The corresponding plot for S xx (q, u>) and qh=\n at T=T C 
is displayed in Fig. 10. The dynamic structure factors 
for the largest three lattice sizes fall onto a single curve 
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FIG. 8: Linewidths (a.)T xx (q) and (b)r zz (q) of the spin-wave 
peak in S kk (q,co) for the [11] direction and L = 300. The 
solid lines are linear fits to the data points. 



within the bounds of their error bars. Our MC-SD data 
clearly show that even above Tkt the dynamic structure 
factor scales quite well, consistent with the results_of the 
2D XY model on a square lattice by Evertz et al&3. For 
the out-of-plane component S zz , we do observe scaling 
(not shown). Thus we can validate the dynamic finite- 
size scaling theory. 

In summary, we have studied the dynamic critical be- 
havior of the classical 2D antiferromagnetic planar(XF) 
model on a triangular lattice using a combination of 
Monte Carlo methods and spin-dynamics techniques 
which use recently developed decomposition methods in 
order to increase the efficiency of the simulations. We 
have calculated the dynamic structure factor S kk (q, w) at 
temperatures below, at, and above T c with lattice sizes 
L < 300. 

At T < Tkt, very strong, sharp spin- wave peaks occur 
in both the in-plane and out-of-plane components of the 
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dynamic structure factor S(q,u>). As T increases, the 
spin-wave peak widens and its position moves towards 
lower to. In addition to the spin-wave peak, S xx (q, to) dis- 
plays a central peak and additional two-spin- wave peaks. 
No central peaks are visible in S zz {q,u>). Above Tkt, 
S xx (q, lo) and S zz (q, lu) apparently display a central peak 
and spin- wave signatures are still in S zz . In a qualitative 
sense, the same behavior was observed in spin-dynamics 
simulations of the 2D XY model on a square latticeE-i 

Below T c , where long-range order appears in the stag- 
gered chirality, we observed an almost dispersionless 
domain- wall peak at high to. The domain- wall peak 
broadens and its intensity decreases as the temperature 
increases. However the intensity of the domain- wall peak 
does not depend on lattice size. The domain-wall peak 
disappears for all q above T c . 



FIG. 9: Finite-size scaling of the dynamic structure factor 
S xx for 2=1.0, qL= 4 tt in the [11] direction, and T=T KT . 
Data points collapse onto the same curve for sufficiently large 
L. Only the data from small lattices deviate systematically. 
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FIG. 10: Finite-size scaling of the dynamic structure factor 
S xx for 2=1.0, <?L= 4 tt in the [11] direction, and T=T C . Data 
points collapse onto the same curve for sufficiently large L. 
Only the data from small lattices deviate systematically. 



The lineshape of the central, spin-wave, and domain- 
wall peaks is captured reasonably well by a Lorentzian 
form. We estimated the dynamic critical exponent z 
from the finite-size scaling of the characteristic frequency 
Lu m . Our determined value is z=l. 002(3) that is in agree- 
ment with-that of the 2D XY model on a square lattice 
z=1.00(4)H 

Finally, we have examined the dynamic finite-size scal- 
ing theory and we found that the characteristic frequency 
cu^(qL) and the dynamic structure factor S kk (q,uj) itself 
scale very well. 
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